/*
Copyright 2004 John Tsiombikas <nuclear@siggraph.org>

This file is part of the graphics core library.

the graphics core library is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.

the graphics core library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with the graphics core library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

/* 3D Curves
 * 
 * author: John Tsiombikas 2003
 * modified:
 * 		John Tsiombikas 2004
 * 		Mihalis Georgoulopoulos 2004
 */

#include <cstdio>
#include <cstring>
#include <cmath>
#include <cctype>
#include <cassert>
#include "curves.hpp"
#include "common/err_msg.h"

Curve::Curve() {
	arc_parametrize = false;
	ease_curve = 0;
	samples = 0;

	set_ease_sample_count(100);
}

Curve::~Curve() {
	delete [] samples;
	
}

void Curve::set_arc_parametrization(bool state) {
	arc_parametrize = state;
}

#define Param	0
#define ArcLen	1

void Curve::sample_arc_lengths() {
	const int samplesPerSegment = 30;
	sample_count = get_segment_count() * samplesPerSegment;

	arc_parametrize = false;	// to be able to interpolate with the original values

	samples = new Vector2[sample_count];
	Vector3 prevpos;
	scalar_t step = 1.0f / (scalar_t)(sample_count-1);
	for(int i=0; i<sample_count; i++) {
		scalar_t t = step * (scalar_t)i;
		Vector3 pos = interpolate(t);
		samples[i][Param] = t;
		if(!i) {
			samples[i][ArcLen] = 0.0f;
		} else {
			samples[i][ArcLen] = (pos - prevpos).length() + samples[i-1][ArcLen];
		}
		prevpos = pos;
	}

	// normalize arc lenghts
	scalar_t maxlen = samples[sample_count-1][ArcLen];
	for(int i=0; i<sample_count; i++) {
		samples[i][ArcLen] /= maxlen;
	}

	arc_parametrize = true;
}

static int binary_search(Vector2 *array, scalar_t key, int begin, int end) {
	int middle = begin + ((end - begin)>>1);

	if(array[middle][ArcLen] == key) return middle;
	if(end == begin) return middle;

	if(key < array[middle][ArcLen]) return binary_search(array, key, begin, middle);
	if(key > array[middle][ArcLen]) return binary_search(array, key, middle+1, end);
	return -1;	// just to make the compiler shut the fuck up
}

scalar_t Curve::parametrize(scalar_t t) const {
	if(!samples) const_cast<Curve*>(this)->sample_arc_lengths();

	int samplepos = binary_search(samples, t, 0, sample_count);
	scalar_t par = samples[samplepos][Param];
	scalar_t len = samples[samplepos][ArcLen];

	// XXX: I can't remember the significance of this condition, I had xsmall_number here
	// previously and if t was 0.9999 it broke. I just changed the number blindly which fixed
	// the breakage but I should investigate further at some point.
	if((len - t) < 0.0005) return par;

	if(len < t) {
		if(!samplepos) return par;
		scalar_t prevlen = samples[samplepos-1][ArcLen];
		scalar_t prevpar = samples[samplepos-1][Param];
		scalar_t p = (t - prevlen) / (len - prevlen);
		return prevpar + (par - prevpar) * p;
	} else {
		if(samplepos >= sample_count) return par;
		scalar_t nextlen = samples[samplepos+1][ArcLen];
		scalar_t nextpar = samples[samplepos+1][Param];
		scalar_t p = (t - len) / (nextlen - len);
		return par + (nextpar - par) * p;
	}

	return par;	// not gonna happen
}


#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define MAX(a, b) ((a) > (b) ? (a) : (b))

scalar_t Curve::ease(scalar_t t) const {
	if(!ease_curve) return t;

	const_cast<Curve*>(this)->ease_curve->set_arc_parametrization(true);
	scalar_t et = ease_curve->interpolate(t).y;

	return MIN(MAX(et, 0.0f), 1.0f);
}


void Curve::add_control_point(const Vector3 &cp) {
	control_points.push_back(cp);
	delete [] samples;
	samples = 0;
}

void Curve::remove_control_point(int index) {
	if(index < 0 || index >= control_points.size()) return;

	ListNode<Vector3> *node = control_points.begin();
	for(int i=0; i<index; i++) {
		node = node->next;
	}

	control_points.erase(node);
}

#define MIN(a, b)	((a) < (b) ? (a) : (b))
#define MAX(a, b)	((a) > (b) ? (a) : (b))

Vector3 *Curve::get_control_point(int index) {
	index = MAX(0, MIN(index, control_points.size() - 1));

	ListNode<Vector3> *node = control_points.begin();
	for(int i=0; i<index; i++) {
		assert(node);
		node = node->next;
	}
	
	assert(node);
	return &node->data;
}

int Curve::get_point_count() const {
	return control_points.size();
}

void Curve::set_ease_curve(Curve *curve) {
	ease_curve = curve;
}

void Curve::set_ease_sample_count(int count) {
	ease_sample_count = count;
	ease_step = (int)(1.0f / (scalar_t)ease_sample_count);
}

Vector3 Curve::operator ()(scalar_t t) const {
	return interpolate(t);
}

///////////////// B-Spline implementation ////////////////////

int BSplineCurve::get_segment_count() const {
	return control_points.size() - 3;
}

Vector3 BSplineCurve::interpolate(scalar_t t) const {
	if(t > 1.0) t = 1.0;
	if(t < 0.0) t = 0.0;

	if(control_points.size() < 4) return Vector3(0, 0, 0);

	if(arc_parametrize) {
		t = ease(parametrize(t));
	}

	// find the appropriate segment of the spline that t lies and calculate the piecewise parameter
	t = (scalar_t)(control_points.size() - 3) * t;
	int seg = (int)t;
	t -= (scalar_t)floor(t);
	if(seg >= get_segment_count()) {
		seg = get_segment_count() - 1;
		t = 1.0f;
	}
	
	ListNode<Vector3> *iter = const_cast<BSplineCurve*>(this)->control_points.begin();
	for(int i=0; i<seg; i++) iter = iter->next;

	Vector3 Cp[4];
	for(int i=0; i<4; i++) {
        Cp[i] = iter->data;
		iter = iter->next;
	}

	Vector3 res;
	res.x = bspline(Cp[0].x, Cp[1].x, Cp[2].x, Cp[3].x, t);
	res.y = bspline(Cp[0].y, Cp[1].y, Cp[2].y, Cp[3].y, t);
	res.z = bspline(Cp[0].z, Cp[1].z, Cp[2].z, Cp[3].z, t);

	return res;
}

//////////////// Catmull-Rom Spline implementation //////////////////

int CatmullRomSplineCurve::get_segment_count() const {
	return control_points.size() - 1;
}

Vector3 CatmullRomSplineCurve::interpolate(scalar_t t) const {
	if(t > 1.0) t = 1.0;
	if(t < 0.0) t = 0.0;

	if(control_points.size() < 2) return Vector3(0, 0, 0);

	if(arc_parametrize) {
		t = ease(parametrize(t));
	}

	// find the appropriate segment of the spline that t lies and calculate the piecewise parameter
	t = (scalar_t)(control_points.size() - 1) * t;
	int seg = (int)t;
	t -= (scalar_t)floor(t);
	if(seg >= get_segment_count()) {
		seg = get_segment_count() - 1;
		t = 1.0f;
	}

	Vector3 Cp[4];
	ListNode<Vector3> *iter = const_cast<CatmullRomSplineCurve*>(this)->control_points.begin();
	for(int i=0; i<seg; i++) iter = iter->next;

	Cp[1] = iter->data;
	Cp[2] = iter->next->data;
	
	if(!seg) {
		Cp[0] = Cp[1];
	} else {
		Cp[0] = iter->prev->data;
	}
	
	if(seg == control_points.size() - 2) {
		Cp[3] = Cp[2];
	} else {
		Cp[3] = iter->next->next->data;
	}

	Vector3 res;
	res.x = catmull_rom_spline(Cp[0].x, Cp[1].x, Cp[2].x, Cp[3].x, t);
	res.y = catmull_rom_spline(Cp[0].y, Cp[1].y, Cp[2].y, Cp[3].y, t);
	res.z = catmull_rom_spline(Cp[0].z, Cp[1].z, Cp[2].z, Cp[3].z, t);

	return res;
}

/* BezierSpline implementation - (MG)
 */
int BezierSpline::get_segment_count() const
{
	if (control_points.size() < 0) return 0;
	return control_points.size() / 4;
}

Vector3 BezierSpline::interpolate(scalar_t t) const
{
	if (!get_segment_count()) return Vector3(0, 0, 0);

	if (arc_parametrize)
	{
		t = ease(parametrize(t));
	}
 
	t = (t * get_segment_count());
	int seg = (int) t;
	t -= (scalar_t) floor(t);
	if (seg >= get_segment_count())
	{
		seg = get_segment_count() - 1;
		t = 1.0f;
	}

	seg *= 4;
	ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
	for (int i = 0; i < seg; i++) iter = iter->next;
	
	Vector3 Cp[4];
	for (int i = 0; i < 4; i++)
	{
		Cp[i] = iter->data;
		iter = iter->next;
	}
	
	// interpolate
	return bezier(Cp[seg], Cp[seg + 1], Cp[seg + 2], Cp[seg + 3], t);
}

Vector3 BezierSpline::get_tangent(scalar_t t)
{	
	if (!get_segment_count()) return Vector3(0, 0, 0);

	if (arc_parametrize)
	{
		t = ease(parametrize(t));
	}
 
	t = (t * get_segment_count());
	int seg = (int) t;
	t -= (scalar_t) floor(t);
	if (seg >= get_segment_count())
	{
		seg = get_segment_count() - 1;
		t = 1.0f;
	}

	seg *= 4;
	ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
	for (int i = 0; i < seg; i++) iter = iter->next;
	
	Vector3 Cp[4];
	for (int i = 0; i < 4; i++)
	{
		Cp[i] = iter->data;
		iter = iter->next;
	}
	
	// interpolate
	return bezier_tangent(Cp[0], Cp[1], Cp[2], Cp[3], t);
}

Vector3 BezierSpline::get_control_point(int i) const
{	
	
	ListNode<Vector3> *iter = const_cast<BezierSpline*>(this)->control_points.begin();
	for (int j = 0; j < i; j++) 
	{
		if (!iter->next)
		{
			return Vector3(0, 0, 0);
		}
		iter = iter->next;
	}
	
	return iter->data;
}



bool save_curve(const char *fname, const Curve *curve) {
	FILE *fp = fopen(fname, "w");
	if(!fp) {
		error("failed to save the curve %s", curve->name.c_str());
		return false;
	}

	fputs("curve_3dengfx\n", fp);
	fputs(curve->name.c_str(), fp);
	fputs("\n", fp);

	if(dynamic_cast<const BSplineCurve*>(curve)) {
		fputs("bspline\n", fp);
	} else if(dynamic_cast<const CatmullRomSplineCurve*>(curve)) {
		fputs("catmullrom\n", fp);
	} else if(dynamic_cast<const BezierSpline*>(curve)) {
		fputs("bezier\n", fp);
	} else {
		error("unknown spline type, save failed %s", curve->name.c_str());
		fclose(fp);
		remove(fname);
		return false;
	}

	fprintf(fp, "%d\n", curve->control_points.size());

	const ListNode<Vector3> *node = curve->control_points.begin();
	while(node) {
		fprintf(fp, "%f %f %f\n", node->data.x, node->data.y, node->data.z);
		node = node->next;
	}

	fclose(fp);
	return true;
}

Curve *load_curve(const char *fname) {
	FILE *fp = fopen(fname, "r");
	if(!fp) {
		error("failed to open file %s", fname);
		return 0;
	}

	char buffer[256];

	fgets(buffer, 256, fp);
	if(strcmp(buffer, "curve_3dengfx\n") != 0) {
		error("load_curve failed, %s is not a curve file", fname);
		fclose(fp);
		return 0;
	}

	Curve *curve;

	fgets(buffer, 256, fp);
	std::string name = buffer;

	fgets(buffer, 256, fp);
	if(!strcmp(buffer, "bspline\n")) {
		curve = new BSplineCurve;
	} else if(!strcmp(buffer, "catmullrom\n")) {
		curve = new CatmullRomSplineCurve;
	} else /*if(!strcmp(buffer, "bezier"))*/ {
		error("unsupported curve type (%s) or not a curve file", buffer);
		fclose(fp);
		return 0;
	}

	curve->name = name;

	fgets(buffer, 256, fp);
	if(!isdigit(buffer[0])) {
		error("load_curve failed, %s is not a valid curve file (count: %s)", fname, buffer);
		delete curve;
		fclose(fp);
		return 0;
	}
	int count = atoi(buffer);

	int failed = count;
	for(int i=0; i<count; i++, failed--) {
		fgets(buffer, 256, fp);
		if(!isdigit(buffer[0]) && buffer[0] != '.' && buffer[0] != '-') {
			break;
		}
		float x = atof(buffer);

		char *ptr = strchr(buffer, ' ');
		if(!ptr || (!isdigit(ptr[1]) && ptr[1] != '.' && ptr[1] != '-')) {
			break;
		}
		float y = atof(++ptr);
		
		ptr = strchr(ptr, ' ');
		if(!ptr || (!isdigit(ptr[1]) && ptr[1] != '.' && ptr[1] != '-')) {
			break;
		}
		float z = atof(++ptr);

		curve->add_control_point(Vector3(x, y, z));
	}

	fclose(fp);

	if(failed) {
		error("load_curve failed to read the data, %s is not a valid curve file", fname);
		delete curve;
		return 0;
	}

	return curve;
}
